Code
from pyrosm import OSM
from pyrosm import get_data
import matplotlib.pyplot as plt
import osmnx as ox
import geopandas as gpd
from osgeo import gdal
import pygeoprocessing as pygeoAdvanced Geocomputation
imports
| Spec | Value (H, M, L) | Effort (H, M, L) |
|---|---|---|
| Get OSM data | H | L |
| Get remote sensing data | H | M |
| Create class correspondence | M | M |
| Rasterize OSM data | H | H |
| Compare OSM and remote sensing rasters | H | M |
Review the spec list. Cross of features that are too high effort for too low value until you have a list of features that can be implemented in 2 weeks. That is your MVP! Success: Now that you have your MVP spec. What is the key metric for success? If you achieve XYZ you know that your MVP was successful. Write 1 sentence outlining your 1 key metric. Reflection: After you implement your MVP, write 6-7 sentences reflecting on your experiences in developing and running your MVP. Was it successful? What did you learn? Would you have done something differently if you were to repeat the project?
from pyrosm import OSM
from pyrosm import get_data
import matplotlib.pyplot as plt
import osmnx as ox
import geopandas as gpd
from osgeo import gdal
import pygeoprocessing as pygeo# Get the geographic boundary for St. Louis County, Minnesota
gdf = ox.geocode_to_gdf("Lakeville, Minnesota, USA")
gdf_path = '/Users/mbraaksma/Files/advgeocomp2024/advgeocomp2024/mvp01/data/lakeville.gpkg'
# albers_crs = {
# 'proj': 'aea', 'lat_1': 29.5, 'lat_2': 45.5, 'lat_0': 23,
# 'lon_0': -96, 'x_0': 0, 'y_0': 0, 'datum': 'WGS84', 'units': 'm', 'no_defs': True
# }
# gdf = gdf.to_crs(albers_crs)
gdf.to_file(gdf_path, driver='GPKG')# File paths
input_raster = '/Volumes/T7 Touch/Research/Reservation_Land/NLCD/raw/nlcd_2021_land_cover_l48_20230630/nlcd_2021_land_cover_l48_20230630.img'
input_vector = '/Users/mbraaksma/Files/advgeocomp2024/advgeocomp2024/mvp01/data/lakeville.gpkg'
output_raster = '/Users/mbraaksma/Files/base_data/lulc/nlcd/nlcd_clipped.tif'
# Open the raster
raster = gdal.Open(input_raster)
# Use gdal.Warp to clip raster with the vector file
gdal.Warp(
output_raster,
raster,
format='GTiff',
cutlineDSName=input_vector, # Path to vector file
cropToCutline=True, # Only crop the raster to the extent of the vector
dstNodata=-9999, # Set nodata value for the output raster
)
# Close the dataset
raster = NoneWarning 1: for band 1, destination nodata value has been clamped to 0, the original value being out of range.
gdf = ox.features_from_place("Lakeville, Minnesota, USA", tags = {"landuse": True})
gdf.explore(column='landuse')/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:
invalid value encountered in intersects
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/set_operations.py:340: RuntimeWarning:
invalid value encountered in union
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:
invalid value encountered in intersects
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:
invalid value encountered in intersects
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:
invalid value encountered in intersects
/Users/mbraaksma/mambaforge/envs/geovenv1/lib/python3.10/site-packages/shapely/predicates.py:798: RuntimeWarning:
invalid value encountered in intersects
lulc_path = '/Users/mbraaksma/Files/base_data/lulc/nlcd/nlcd_clipped.tif'
lulc_data = gdal.Open(lulc_path)
lulc_array = lulc_data.GetRasterBand(1).ReadAsArray()
# Close the dataset
lulc_data = None
# Plot the raster data
plt.figure(figsize=(10, 8))
plt.imshow(lulc_array, cmap='viridis') # Adjust colormap if needed
plt.colorbar(label='Raster Values')
plt.title('Raster Plot')
plt.show()